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The Biham-Middleton-Levine (BML) traffic model, a cellular automaton with east-bound and 
north-bound cars moving by turns on a square lattice, has been an underpinning model in the study 
of collective behaviour by cars, pedestrians and even internet packages. Contrary to initial beliefs 
that the model exhibits a sharp phase transition from freely flowing to fully jammed, it has been 
reported that it shows intermediate stable phases, where jams and freely flowing traffic coexist, 
but there is no clear understanding of their origin. Here, we analyze the model as an anisotropic 
system with a preferred fluid direction (north-east) and hnd that it exhibits two differentiated phase 
transitions: either if the system is longer in the flow direction (longitudinal) or perpendicular to 
it (transversal). The critical densities where these transitions occur enclose the density interval of 
intermediate states and can be approximated by mean-held analysis, all derived from the anisotropic 
exponent relating the longitudinal and transversal correlation lengths. Thus, we arrive to the in¬ 
teresting result that the puzzling intermediate states in the original model are just a superposition 
of these two different behaviours of the phase transition, solving by the way most mysteries behind 
the BML model, which turns to be a paradigmatic example of such anisotropic critical systems. 
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In the recent urbanization era, society faces an in¬ 
evitably increase of traffic congestion, turning its at¬ 
tention onto urban road networks. Nowadays is widely 
assumed that a proper understanding of the mecha¬ 
nisms leading jamming processes is indispensable for im¬ 
proving the efficiency of transportation systems. The 
Biham-Middleton-Levine (BML) model [1] is, perhaps, 
the simplest traffic cellular automaton able to exhibit 
self-organization, pattern formation and phase transi¬ 
tions m- Hence, much extensive research on hux an 
collective behavior has been based on it, not just for car 
traffic (see ISIE] and Refs, in di) but also for pedes¬ 
trian 0 uni and information packages on the Internet 
m- For more than a decade, it has been believed that 
at a certain critical car density pc, the system exhibits 
what seems to be a hrst order phase transition between 
two phases: a free-howing phase, where all cars move 
freely at all time steps (the average velocity of cars v=l) 
and a completely jamming phase, where no car moves at 
all (u=0). The value of pc decreases with increasing sys¬ 
tem size, possibly reaching the value Pc=0 as the system 
size approaches infinity. Then, it was thought that the 
BML model would be similar to other well known sys¬ 
tems in statistical physics exhibiting phase transitions, 
e.g. percolation. However, all these conventional be¬ 
liefs started being reconsidered, since W.K. Young [12] 
and R. M. D’Souza [131 IHl realized the existence of in¬ 
termediate stable phases where free-flowing and jamming 
phases coexist (see FigH- The structure of these states 
is highly regular, with jams’ wavefronts moving through 
freely flowing traffic in a wide density region, and the 


value of the average velocity (0<u< 1) is extremely sen¬ 
sitive to the aspect ratio of the underlying lattice. Thus, 
instead of a phase transition as a function of car density, 
the system would exhibit two bifurcation points, limit¬ 
ing a region where intermediate phases would coexist be¬ 
tween the two conventional phases above. Even though 
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FIG. 1. (color online) (a) Average velocity {v) vs density p 
(circles) for a BML model on a 45^-rotated 256 x 256 lattice 
(as depicted in (a), with east-bound cars (brown) and north¬ 
bound cars (black)). The mean-held approach for longitudinal 
systems (red line) and the numerical prediction for an inhnite 
transverse system (blue line) are also included for comparison. 
Typical conhgurations for free how (c), intermediate states 
(d), one global jam (e) and random jams (f) are also included. 
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these states have been described and their asymptotic 
speeds have been predicted m, the exact locations of 
the bifurcation points are very difficult to determine and 
no one knows what truly happens as the system size goes 
to infinity. Moreover, the origin of the intermediate states 
remain as a unsolved puzzle. 

In this Letter we focus on the clear existence of a 
preferred direction in the model dynamics. Most pre¬ 
vious studies have overlooked this feature, despite - as 
we will show - it is the key for unraveling the puzzling 
intermediate states of the BML model. The presence 
of anisotropy has been fundamental in the analysis of 
force networks in granular matter [151 IS] flocking 
in collective animal behavior m, and classical theo¬ 
retical models with this features are the next-nearest- 
neighbor Ising model (ANNNI)[T8] and the driven lat¬ 
tice gas model [19], and it should be naturally reflected 
in the phase transition of BML, making it in principle 
more akin to the anisotropic equivalent in percolation: 
directed percolation (DP) [20l[2T]. By performing a di¬ 
rect anisotropic scaling analysis on the BML model phase 
transition [22l [23] . we found that the jamming process 
behaves distinctively as two separated phase transitions 
along different directions, namely, whether the system is 
longer in the direction of traffic flow (longitudinal sys¬ 
tem) or in the transversal direction (transversal system). 
Our main result is that the puzzling intermediate states 
in the BML model on square lattices emerge just as a 
superposition of these two different transitions, solving - 
by the way - most mysteries behind the BML model. 

Model - The original BML model considers two types 
of cars: east-bound (brown) and north-bound (black), 
moving on a two-dimensional square lattice with peri¬ 
odic boundary conditions. Each lattice site is in one of 
three states: empty, occupied by a brown car, or occu¬ 
pied by a black one. The cars are initially randomly 
distributed over the lattice sites with spatial density p 
(usually taken to be the same, p/2, for both north- and 
eastbound cars). The fully deterministic dynamics is as 
follows: On even (odd) steps, all eastbound (northbound) 
cars synchronously attempt to advance one lattice site to¬ 
ward the east (north). If the site eastward (northward) 
of a car is currently empty, it advances. Otherwise, it re¬ 
mains stationary. The system exhibits, therefore, a pre¬ 
ferred north-east direction for the flux. 

To study the system as an anisotropic one, we ex¬ 
plicitly rotated 45 degrees the lattice, so we could con¬ 
trol the system lengths along the longitudinal (Ly) and 
transversal {L±) directions to the car flow. Although 
this rotation changes the boundary conditions, the sys¬ 
tem still exhibits the same three phases observed in the 
original BML model (see FigQ. In anisotropic sys¬ 
tems, clusters show different correlation lengths along 
the longitudinal and transversal directions, and ^_l, 
respectively, which scale with different exponents as 
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FIG. 2. (color online). Longitudinal and transversal 
correlation lengths from final configurations at densities p in 
the range [0.49-0.54] for square lattices of different sizes. The 
power-law fit gives ^=1.64(3) for the anisotropy exponent. 


i\\'^{p-pc) "II and C-L~(p-Pc) "-L [inilll]- An 
anisotropy exponent, relating the different scaling of the 
two correlation lengths, is defined as the ratio 0=^. Ac¬ 
cording to |22ll24] . when the longitudinal and transver- 
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sal lengths are related by the constraint L_\_ ^ I/|p the 
system behaves as it were effectively isotropic, and stan¬ 
dard finite-size scaling (FSS) theory applies again for all 
percolation quantities, just in terms of the length scale 
Ly. Especially the transition width and the percolation 
threshold scale like 





and 



{pciLiU^)) 



( 1 ) 

The exponent 0 can be estimated numerically from the 
fact that, close to the critical point, the two correlations 
lengths must be related by ^y ^ The symbol () de¬ 
notes averages over final jammed configurations starting 
from different random initial conditions. With this the¬ 
ory in mind, let us define the parallel (perpendicular) 
spatial correlation function [2] as 


G\\(±){P) = , 


( 2 ) 


where cr(x)=l(0) if the site with position x is occu- 
pied(empty), N is the total number of cars and r' is a 
vector in the direction || (T) you want to compute the 
correlation function along. The correlation functions are 
fitted with exponentials Gy(^) oc exp(—r/^y(^)) to esti¬ 
mate C||(±)- Figure presents the correlations lengths 
computed from final configurations of the original BML 
model for different square lattice sizes at densities close 
to the threshold transition, averaging over 50 configura¬ 
tions for each point. A linear regression yields an esti¬ 
mate of the anisotropy exponent ^=1.64 ±0.03 (Here and 
everywhere the error bars are la). 

With the exponent 0 in hand, we studied the phase 

transition for both longitudinal and transversal sys- 

1 

terns. Longitudinal ones ran on lattices with sizes x 
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FIG. 3. (color online). Finite-size scaling analysis for the phase transition of both longitudinal and transversal systems, (a) 
Transition curves for several longitudinal systems, including the mean-field prediction (magenta) and a typical configuration 
for the jammed phase (inset), (b) Transition curves for several transversal systems (sizes as in (a), but with Ly and L± 
interchanged), (c) Scaling of A“^(L||) for both cases, giving z/y =3.0(3) and u±=1.8{2). (d) Scaling of the finite critical densities 
give pc|| =0.521(2) and =0.283(2). All systems were simulated for at least 1.5 x 10® time steps or until convergence (r'^O or 


v—1). 


I/||. Figure ^ a) shows that velocity begins to decrease 
smoothly with increasing p until the abrupt onset of full 
jamming (n=0) at a certain pcy • Right there, an approx¬ 
imately uniform distribution of jams spans through the 
whole system (Figj^a), inset). With Eq. 0 in mind, the 
FSS analysis (FigMc) and (d)) gives pcy =0-521(2) and 
z/|| =3.0(3), which would imply v^=1.8{2). 

Quite differently, transversal systems (i.e., lattices with 
1 

L\\ X I/|p in the gridlock phase show a single and well 
localized jam on an empty background (FigJ^b)). These 
systems exhibit a sharp phase transition between free¬ 
flow and a completely gridlock phases, but at a lower 
density pcj_- The FSS analysis (also in FigJ^c) and (d)) 
gives pc^ =0.283(1) and =2.2(1), in agreement with 
the previous result. Hence, the intermediate states in 
square lattices [13] emerge just as a consequence of the 
combination of these two phase transitions (FigJ^a) and 

(b)). 

The critical densities for both longitudinal and traverse 
transitions can be approximated by using a mean-field 
analysis, inspired by [25]. Consider the mean velocity 
of brown cars (by symmetry, the reasoning is also valid 
for black cars). A brown car will stop either because 
it is crossed by a black car or because it queues behind 
another brown car. At a random initial configuration, 
the probability that a car is crossed or queued is p^, that 
is, at the beginning of the simulation the proportion of 
stopped cars Pstop must be equal to p. Let us define 
(c^^) as the proportion of stopped cars that are crossed 
(queued). In previous works [25] [26], it has been assume 


that c^^=c^^=0.5, but we will see that this is not the 
case. If Pstop ^ P for some time steps, the probability a 
cell to be occupied by a stopped crossed (queued) car will 
be Since black cars spend on average a 

time 1/v on a site, they will reduce the speed of brown 
cars from unity by jv. Similarly, the extra amount 

of time that a brown car stays on a site will be given 
by ^ — 1, reducing the average speed by c^^p^(^ ~ 1)- 
Hence, a self-consistency equation for the average speed 
V will be 



which gives pc as the critical density at which the equa¬ 
tion ceases to give a real solution. 

Consider longitudinal systems first, where a uniform 
distribution of longitudinal jams arise. Before solving 
Eq. 0, let us study the evolution on time of the quan¬ 
tities involved. On one hand, Eigj^a) shows the time 
evolution of Pstop for different values of the car density 
for a longitudinal system. At start, Pstop=P for every den¬ 
sity, as expected. Later on, Pstop declines (grows) for low 
(high) densities, and the crossover occurs near the critical 
density pc. Actually, Pstop after one time step equals p for 
p ^ 0.57 in a large system (Eigj^b)). On the other hand, 
c^>^=c^^=0.5 only at start, but after just one time 
step they change to 600(5) and c^^^O.400(5), 

remaining there for some time steps (EigQc)). By re¬ 
placing these two values into Eq. 0, the mean-field ap¬ 
proach predicts a critical density py =0.563(5), in agree- 
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FIG. 4. (color online). Time evolution of the fraction of 
stopped cars Pstop (a) and of the fraction of crossed (c^|) and 
queued (c^^) stopped cars (b) in a 1438 x 82 longitudinal 
system. The value of Pstop after one time step for several 
densities is also shown (c). (d) Schematic representation of 
the jamming process for transversal systems 


ment with the value obtained from both finite-size scaling 
and the crossover described above. 

Transversal systems, in contrast, collapse into a single 
big jam. Thus, the previous reasoning is impractical here, 
because the local density becomes inhomogeneous during 
the collapse: increasing in the vicinity of the growing jam 
and dropping to zero elsewhere. Let us consider that the 
density at start is low enough to assume that all cars 
are moving. At a certain evolution time the cars have 
condensed into a thinner region (red area in Figj^d)) of 
area Ared ? where the density has reached the critical value 
PcII for longitudinal systems. All cars in this region come 
from a larger region of area Agreen in the initial configura¬ 
tion (green area in FigQd)). Let us assume that the red 
area is a fraction of the green one equals to the fraction of 
cars that have stopped, AredMgreen=Pstop- Because the 
number of cars does not change, at the critical density 

_ ^red _ 2 f a\ 

Pgreen Pred~Z Pcy ’ 

^green 

where we have used the result that Pstop = p at the crit¬ 
ical density for longitudinal systems. Then, the critical 
density for transversal systems will be PcjL—p‘c\\ =0.271(3), 
which approximates quite well the value pcj^ =0.283(2) 
from the simulations (FigJ^d)). 

The values for and can also be estimated 

from the scaling analysis. Consider a typical jam (as the 
one in FigJ^a)). It consists of two sets of queues, east- 
bound (brown) and north-bound (black), intersecting 
along a longitudinal curve. Crossed cars localize on this 
intersection, and its number is proportional to 
Queued cars fill the whole jam, and its number is propor¬ 
tional to the jam’s area (i.e. to xLj_). Thus, 

and . Replacing into + c^^=l. 


they give c^^=0.582(1) and c^^=0.418(1), very close to 
the simulation results. Using these values in the mean- 
field equation gives py =0.5672(2), also close to the sim¬ 
ulations. 

Taken together, our results suggest that the puzzling 
intermediate states in the BML model on square lattices 
are, actually, a consequence of the system’s anisotropy, 
which produces two different phase transitions: one for 
transverse systems, with pc^ =0.283(1), and another for 
longitudinal ones, with pc^ =0.521(1). Indeed, the first 
critical density corresponds with the lower bifurcation 
point on square lattices, reported at p=0.315 m , con¬ 
tradicting the general believe that this would go to zero 
for infinite systems; similarly, the second critical den¬ 
sity perfectly matches with the value of p=0.52 reported 
as the transition point between self-organized jams and 
random jams [ 211 ] within the conventional understand¬ 
ing of the original BML model. The asymptotic limit of 
the curves from simulations on transverse systems (blue 
line in FigQ and the mean-field curve for longitudinal 
systems (Eq. (§, red line in FigQ perfectly enclose the 
zone of intermediate states. 

Despite having studied the critical features of the BML 
phase transitions with the DP formalism, the obtained 
critical exponents are not compatible with the universal¬ 
ity class of the directed percolation. Instead, the critical 
behavior here coincides with the reported for the parity- 
conserving universality class (pc)iiniiii]. The origin of 
this relationship remains as an open question for future 
research. 

By finding the origin of the intermediate states of the 
BML model, we have built a very complete description of 
such a fundamental model for traffic flow, illustrating at 
the same time the power of the finite-size scaling analysis 
for anisotropic systems, where the BML model seems to 
be a paradigmatic example. We look forward to see many 
more results of this enlightening analysis technique in the 
future. 
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